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The classical problem of the Josephson junction of arbitrary length W in the presence of externally 
applied magnetic fields [H) and transport currents (J) is reconsidered from the point of view of 
stability theory. In particular, we derive the complete infinite set of exact analytical solutions for the 
phase difference that describe the current-carrying states of the junction with arbitrary W and an 
arbitrary mode of the injection of J. These solutions are parameterized by two natural parameters: 
the constants of integration. The boundaries of their stability regions in the parametric plane are 
determined by a corresponding infinite set of exact functional equations. Being mapped to the 
physical plane {H,J), these boundaries yield the dependence of the critical transport current Jc 
on H. Contrary to a wide-spread belief, the exact analytical dependence = Jc (H) proves to be 
multivalued even for arbitrarily small W. What is more, the exact solution reveals the existence of 
unquantized Josephson vortices carrying fractional flux and located near one of the junction edges, 
provided that J is sufficiently close to Jc for certain finite values of H. This conclusion (as well as 
other exact analytical results) is illustrated by a graphical analysis of typical cases. 

PACS numbers; 74.50.-|-r, 03.75.Lm, 02.30.Oz 



Based on mathematical methods of stability theory, we reconsider the classical physical problem^ "^■'^ of current- 
carrying states of the Josephson junction of arbitrary length W in external magnetic fields. Although the problem 
was first posed over four decades ago^i^i^ and ever since has found numerous practical applicationSj^'^'i its complete 
analytical solution has not been obtained in the previous literature. Here, we derive this solution and show that 
it leads to new and important physical conclusions: the multivaluedness of the exact analytical dependence of the 
critical transport current on the applied field for arbitrarily small W, and the existence of unquantized Josephson 
vortices carrying fractional flux. This paper can be considered as a logical continuation of the investigation initiated 
in our preceding publication^^ where we have derived the complete analytical solution for the Josephson junction in 
external magnetic fields in the absence of transport currents. 

To remind the reader of the standard formulation of the problem, we consider the geometry presented in Fig. [T] 
Here, the x axis is perpendicular to the insulating layer I (the barrier) between two identical superconductors S; the 
y axis is along the barrier whose length is W = 2L G (0, 00). A constant, homogeneous external magnetic field H is 
applied along the axis z: H = {0,0, H > 0). Full homogeneity along the z axis is assumed. The transport current J 
is injected along the axis x: J = (J, 0, 0) . 

In the region of field penetration, the electrodynamics of the junction in equilibrium is fully described by a time- 
independent phase difference at the barrier, cj) = (t){y). Using the dimensionless units introduced in Ref4, we can 
write down the local magnetic field and the Josephson current density as^ 
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FIG. 1: The geometry of the problem: t is the thickness of the barrier; W = 2L is the length of the barrier; A is the London 
penetration depth; d = 2A + 1 is the width of the field-penetration region (shaded). The external magnetic field H is directed 
into the plane of the figure, and the transport current J is along the axis x. 



If it is symmetric with respect to the plane (y, z), we have: 
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Solutions to ([3]), (|3]) are supposed to satisfy an obvious physical requirement: they must be stable with respect to any 
infinitesimal perturbations. (Unstable solutions that do not meet this requirement are physically unobservable and 
should be rejected.) 

Unfortunately, the standard boundary- value problem (H]) is mathematically ill-posed;^ (i) for | J| larger than 
certain Jmax — Jmax (H, L) , it does not admit any solutions at all; (ii) aside from stable (physical) solutions, there may 
exist unstable (unphysical) solutions for the same H and J; (iii) for the same H and J, there may exist several different 
physical solutions. An immediate consequence of this ill-posedness is as follows: although the general solution to ([3|) 
is well-known^i^ the constants of integration specifying particular physical solutions cannot be determined directly 
from the boundary conditions (|4|). 

In view of the above-mentioned mathematical difficulties, the previous analysis of the problem ([3]), ^ was concen- 
trated mainly on finding the dependence Jmax = Jmax (H) (for particular values of L) without trying to establish the 
exact analytical form of current-carrying solutions. (It should be noted that the quantity Jmax itself was identified 
with the experimentally observable critical current Jc, i.e., the identity Jmax = Jc was assumed.) 

For the case L ^ 1, there existed*^ a simple analytical approximation for the dependence Jmax = Jmax (Jf) (the 
so-calledii^ "Fraunhofer pattern"). As to the case L > 1, only particular numerical results were obtained. Thus, 
Owen and Scalapinc^ established the dependence Jmax = Jmax (H) only for L = 5: it proved to be multivalued. The 
numerical method of Ref.^ was later employed to study the effect of asymmetric injection of the transport currenti^i 
Unfortunately, all these numerical results could tell very little about general properties of the current-carrying states 
for arbitrary L G (0, oo). Besides, no analytical expressions were derived that could serve for direct determination of 

Jmax- 

On the other hand, attempts were mad o^^'^'^ to simplify the computational procedure by transforming the 
boundary- value problem ([3]), ^ into an equivalent initial- value problem. Although these attempts did not pro- 
duce exact analytical solutions, we note that Refsi^^ii^ introduced a new, more satisfactory mathematical definition of 
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the observable critical current J^- it was identified with the boundary of the stability regions of the current-carrying 
configurations. The same mathematical definition of Jc was employed in Refs.^*'^'^ concerned with certain nontrivial 
generalizations of the boundary- value problem ([3]), (jH). Unfortunately, exact analytical expressions for the physical 
solutions to ([3]), ([U were not found in Refsiiiii^, either. 

As already mentioned, in ReL- we have derived the complete infinite set of exact physical solutions to ([3]), ([4]) under 
the condition J = 0. The approach of Refi^. consists in a certain generalization of the boundary conditions and an 
application of methods of stability theory at an early stage of the consideration. The same approach is adopted in 
this paper for the general case J ^ 0. Thus, we derive a complete set of exact particular solutions to ([3]) that are 
stable under the condition that ^ is fixed at the boundaries y = ±L [for arbitrary L e (0,cxd)]. These solutions 
are parameterized by two natural parameters: the constants of integration of ([3]). The boundaries of their stability 
regions are determined by a corresponding infinite set of exact functional equations. The physical interpretation of 
the obtained solutions stems from the fact that the boundary conditions in the form ([S]), ^ (or their modification 
for the case of asymmetric injection of J) realize a mapping of the stability regions from the parametric plane to the 
physical plane {H, J) . 

In Sec. II, we present a static method of the analysis of stability based on the minimization of the generating free- 
energy functional. A Sturm-Liouville eigenvalue problem that plays a key role in the analysis of stability is discussed. 
In Sec. Ill, we derive the complete set of exact stable analytical solutions to ([3]), ([4]) under the condition H > 0, J > 0. 
A numerical analysis of several typical cases is carried out. In Sec. IV, we elaborate on major physical implications 
of the exact analytical solutions. Graphic illustrations are presented. Generalizations to the case of arbitrary sign of 
H and J, and to the case of asymmetric injection of J are considered. Finally, in Sec. V, we summarize the obtained 
physical and mathematical results and make several concluding remarks. 

In Appendix A, an alternative (dynamic) method of the analysis of stability is presented. In Appendix B, functional 
equations for the stability regions are derived. In Appendix C, a certain special solution of the Sturm-Liouville 
eigenvalue problem is considered. 



II. ANALYSIS OF STABILITY 



The stability of the solutions to ©-([S]) can be analyzed by means of two different methods: a statio^ one, and 
a dynamioi^ one. Although they are fully equivalent mathematically, the static method seems to be more natural 
physically: we therefore discuss it in this section. (For the sake of completeness, we outline the dynamic method in 
Appendix A.) 



A. Minimization of the Gibbs free-energy functional 



The static method is based on the minimization of the generating Gibbs free-energy functional. For the boundary- 
value problem ((31), Q, the corresponding functional (in terms of the dimensionless units, ^ and per unit length along 
the z axis) has the following form: 
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As can be easily seen, the stationarity condition of ([7]), 
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dy 



H,J 



0, 



yields the equation for the phase difference ([3]) and the boundary conditions ([4]). 

Note that the functional ^ with J = is analyzed in Ref.*. Basic properties of functionals of the type ([7|) are also 
discussed in RefsJ^ii^: in particular, all the stationary points of ([7]) are either local minima or saddle pointsJ^ 

In full analogy with the case J = Or- the type of a stationary point (f) — (f){y) obeying Q, (|4]) is determined by the 
sign of the lowest eigenvalue ^ = /Xq of the Sturm-Liouville problem 
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Namely, if /zq < 0, the solution (j) = (j){y) corresponds to a saddle point of ^ {S^g ^ 0). Solutions of this type are 
absolutely unstable and hence unphysical. 

On the contrary, the stable physical solutions (f) — 4){y) that minimize ^ are characterized by /xq > {SVIq > 0). 
The boundaries of the stability regions for these solutions {SflQ > 0) are determined by the condition 

A*o = 0, 

or, equivalently, by the solution ipo = ipo iu) to the boundary-value problem 



dy2 



+ cos(/)(zj)?/;o = 0, ye{-L,L), (10) 

^o{y)^0, ye[-L,L]. (12) 

Equation ([H]) can be transformed into Lame's equation. In certain limiting cases, the eigenvalue /i = /xq (and 
the corresponding eigcnfunction = rpo) of the problem ([8]), ([9]) can be found explicitly by perturbation methods: 
see Appendix C. However, since we will mostly need information about the boundaries of the stability regions, the 
consideration of the main part of this paper is based on the fact that the linear boundary- value problem pTH) - (fT^ is 
exactly solvable. The relevant exact analytical solutions are derived in Appendix B. 

III. CURRENT-CARRYING STATES 

As is well-known, the general solution to ([3]) can be easily obtained using the first integral, 

cos(/) = C, -l<C<oo, (13) 



1 " ' 



2 



dy 



where C is the constant of integration. In Ref.^, we have written down the general solution to ^ in the form 
convenient for applications with the boundary conditions In that paper, solutions parameterized by C S [—1, 1) 
and C e (1, +oo) have been termed solutions of type I and type II, respectively. 

As we have shown ior H ^ 0, J = 0,— all the solutions of type I are absolutely unstable. On the contrary, the 
solutions of type II contain, for iJ ^ 0, J = 0, a subclass of stable solutions. 

The case H ^ 0, J ^ is quite different, because both the classes of solutions (of type I and type II) contain 
subclasses of stable current-carrying solutions. [For example, for J = 2H, we have C — cos (/){—£) < 1, since 
(j> {—L) 7^ mod 27r.] In view of continuous dependence of the left-hand side of (fT3|) on C, stable current-carrying 
solutions of type I in the limit C — s- 1 — should coincide with stable current-carrying solutions of type II obtained 
by the limiting procedure C — > 1 + 0. 

Note that, in what follows, we will employ instead of C a standard parametrization constant fc.* Namely, 



for the solutions of type I, and 



fc2 = i±^, 0<fc<l (14) 



0<fc<l (15) 



1 + C" 



for the solutions of type II. Moreover, in this section, we restrict ourselves to symmetric injection of J [conditions 
(|4])], and to the case H > 0, J > 0. (These restrictions will be removed in Sec. IV.) 

According to the scheme outlined in the Introduction, we start with finding all the solutions to ^ that are stable 
under the condition that ^ is fixed at the boundaries y — zLL. These solutions are parameterized by k and the 
second (additive) constant of integration denoted as f3 (for solution of type I) or a (for solutions of type II). The 
boundaries of the stability regions are determined from the solution to the linear boundary- value problem p0|) - p2| . 
Finally, relations ([5]), ^ are employed to map the stability regions from the parametric planes {k,(3) and (fc, a) to 
the physical plane [H, J) . 
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A. Solutions of type I 

The general form of the sohitions of type I is given by^ 

4'± (y) = (2?^ + 1) ± 2 arcsin [ksn{y — yQ,k)] , n = 0, ±1, . . . , (16) 

where sn u is the Jacobian elliptic sine.-^i The constant of integration yo is subject to the restriction 

- K (k) <yo<K{k), (17) 

with K (fc) being the complete elliptic integral of the first kind,— and the constant of integration k is defined by (I14p . 
Taking into account that A: = in (fT6|) corresponds to absolutely unstable solutions with H = J — 0^- we impose the 
condition 

< fc < 1. (18) 

Mathematically, it is convenient to begin the consideration of the current-carrying solutions of type I with the case 
H = 0, J > 0. The solutions for the case H > 0, J > will be obtained from the solutions for if = 0, J > by the 
introduction of a new parameter. 



1. The case H = 0, J > 
The generalized form of the boundary conditions Q for = 0, J > is given by the relations 

^W--l^(-i^), (19) 
dy dy 

^ (L) = const > 0. (20) 
dy 

Using HI]), we find that j/o = -K (fc) in (dH), whereas (gOl) yields (t>=(t>- [L < 2K (fc)]. Finally, setting n = in (fTB|) . 
we obtain 



(t/) = 2 arccos 



cn {y, fc) 
dn(y,fc) 



(21) 



where cn u and dn u are the Jacobian elliptic cosine and the delta amplitude, respectivelyi^i 
This solution is symmetric with respect to reflection: 

i^y) = (y) . (22) 

It is stable only for 

fc e [fc„ 1] , (23) 

where, according to the results of Appendix B, the boundary of the stability region fcc = fcc (L) is implicitly determined 
by the functional equation 

cn {L, fc,) [^E {L, fc,) + (1 - kl) L]+{1- kfj sn (i, fc,) dn [L, kc) = 0, (24) 

with E {u, fc) being the incomplete elliptic integral of the first kind.^^ [We include the point fc = 1 in the definition of 
the stability region because lim^^i^s = 0, which is an absolutely stable solution for the case H = J — 0.] 
Equation ([M]) can be solved analytically in two limiting cases. In particular, for L <§; 1, the solution is 

fce«^. (25) 

For L ^ 1, equation becomes 

K{kc)~L, (26) 



6 




and the solution is 



FIG. 2: The dependence kc = kc (L) (solid line). The stability region is shaded. 



fee « 1 - 8exp(-2L) . 



For arbitrary L € (0,oo), we present the numerical solution to p4)) in Fig. [2] 

Substituting ([2T|) into ([6|), we arrive at the expression for the current J = J {L, k): 

dn [L, k) 

Note that for L = ^ <C 1 expression ([^ reduces to the expected resultii^i^ii^ 

W 

J « y sin (j)s (0) , 



(27) 



(28) 



(29) 



where, by ^ and ((25l) . 

05 (0) = 2 arccos fee 0, ■ 
According to ([28|l . the dependence Jc = (i) is given by 

rsn (L,A:c) 



dn (i, /cc) 



Thus, for L = -f > 1, we get, using l[26 | . (p7 | . 



Jc«2[l-8exp(-W^)] 



(30) 



(31) 



For arbitrary L € (0, oo), the dependence Jc = Jc{L) is presented in Fig. [3] Although Fig. [3] reproduces the old 
results^ obtained by numerical maximization of J, we want to emphasize a substantial methodological difference: 
the curve Jc ~ Jc (L) in Fig- El is nothing but a mapping by means of (j30p of the boundary of the stability region 
kc = kc [L) in Fig. m 



2. The case H > 0, J > 



For fl" > 0, J > 0, instead of IHl), we have 



d(j> 



(32) 
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FIG. 3: The dependence Jc — Jc (L) for H — (solid line). The stability region is shaded. 



Boundary conditions ([5^ break the symmetry (P^ . Taking into account that in the hmit H ^ we must get (|2ip . 
conditions (PD|) and ([5^ can be satisfied by 



(y) = 2 arccos 



cn {y + /3, fc) 
dn {y + (3, k) 



(33) 



where kc is determined by ([24|) . and /3c G [0, X (fc)). The boundary of the stabihty region (3c — Pc{k) is determined 
(see Appendix B) by the solution to the functional equation 



. (L + /3e, fc) cn {L - Pc, k) [-E [L + P^, k) - E [L - (3^, k) + {l - fc') L] 
+ (1 - fc') [sn {L + (3a, k) cn (L - (3^, fc) dn (L + (3c, fc) 
+sn [L - /3c, fc) cn (L + /3c, fc) dn (L - /3c, fc)] =0, fc € [fcc, 1) , 



(34) 



under the condition /3c (fcc) = 0. 

In Fig. m we present the stability region of obtained by numerical evaluation of Eq. for several different 
values of L: L = 0.3 (a "small" junction), L = 1 (a "medium" junction), and L = 3 (a "large" junction). As we can 
see, limfc^i (3c (fc) oo. The asymptotics of /3c (fc) for fc ^ 1 can be established analytically. 

Let us make the substitution 



(3c (fc) - 7c 



(35) 



in Eq. p4p . By proceeding to the limit fc = 1, we obtain a functional equation that determines the dependence 
7c = 7c (L) for fc = 1: 

1 2 

L sinh (L — 7c) sinh {L + jc) — sinh (L — 7c) sinh (L + 7c) cosh {L — 7c) 

1 2 

— -sinh {L + 7c) sinh {L — 7c) cosh (i + 7c) — sinh {L + 7c) cosh {L — 7c) — sinh (L — 7c) cosh {L + 7c) — 0. (36) 

The numerical solution to this equation is given in Fig. \5[ [Note that 7c (L) « L for L ^ 1.] Taking into account 
relation l[35|) . we arrive at the sought asymptotics of /3c (fc) for fc — > 1: 



1 16 
(3cik)»-\nj-^^jciL). 

Accordingly, the limiting form of the current-carrying solution p3p is 

lim (j>s (y) = ((ii {y) = 4 arctan [exp (y - 7)] , 



(37) 



(38) 



where 7 e [7c, 00) (see Fig. [S]). 
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a, (3 
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FIG. 4: The stability regions of and (pp {p — 0,1,2 . . .) in the parametric plane (shaded) for L = 0.3, 1, 3. The dependencies 
l3c = l3c (fc) and Qc = Oc (k) are given by the dashed line and the solid lines, respectively. 



Equations (O and upon the substitution of ([55]) . yield 



fc2 



sn(i + /3. A:) sn {L~P,k) 



dn (L + /3, k) 
sn{L + /3,k) 



dn (L - P, k) 
sn (i — (3, k) 



dn (L + (3, k) dn (L - /3, /c) 
fee [feci), /3G[0,/3e] 



(39) 
(40) 
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L 

FIG. 5: The dependence 7c = 7c (L) (solid line). The stability region is shaded. 



In the limit fc = 1, these relations take the form 



H 



cosh L cosh 7 



J 



cosh (7 — L) cosh 7 ' 
2 sinh L sinh 7 

7 £ [7c, 00) . 



cosh (7 — i) cosh 7 ' 



(41) 
(42) 



By setting f] = (3c in ([M]) and (HDl) , we can obtain a relevant part of the dependence Jc = Jc (H) for arbitrary 
L G (0,00): see Sec. IV. 

We conclude the discussion of the solution 0s by presenting an explicit analytical expression for the special case 
J = 2H that was the subject of numerical evaluation in Refii^. From and (PO]). we find: (3 — L. Substitution 
into (|33p immediately yields 



4>s {y)\j-_ 



:2H 



— 2 arccos 



cn (y + L, k) 



dn (y + L, k) 



k G \krai 1] ; 



(43) 



where km is determined by the condition /3c (fcm) = 



B. Solutions of type II 

We start with the stable type-II solutions for the case H > 0, J — Qi^ ii^iis 

,/,p(y) = 7r(p- 1) + 2am(^| +X(/c) ,fc) , p ^ 2m (rn = 0,1,...) 

(f)p{y) — TTP + 2am (^^,k^ , p = 2m + 1 (to = 0, 1, . . .) 

where am u is the Jacobian amplitude. The stability regions of (|44p . (|45)) are given by 

p^O : k e (fci, 1) ; p = 1,2, . . . : fc e (fcp+i, fcp] , 
where the points k = kp (p = 1, 2, . . .) are the roots of the equations 

pkpK (kp) ^ L, p=l,2, 



(44) 
(45) 

(46) 
(47) 



Solutions (|44| . (|45|) form an infinite set, and the union of their stability regions (|46)) (they interchange for even and 
odd p) is equal to the whole fc-interval (0, 1). The meaning of the parameter p = 0, 1, 2, ... is revealed by the relation 



P = 



(j)p (L) - cj)p (-L) 
27r 



(48) 
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where [. . .J stands for the integer part of the argument. Note also the symmetry property: 



(49) 



For H > 0, J > 0, current-carrying type-II solutions obey the generalized boundary conditions (PO)) . ([5^ that break 
the symmetry (|49|) . These conditions can be satisfied if, in ([44]) and ([45]), we make a shift of the argument y y + ka 
(with a being a new parameter): 



(l)piy) = 7r{p- 1) + 2a,m(Jj- + K (k) +a,kj , a S [0, aj , p = 2m (m = 0, 1, . . .) ; 



(t>p{y) — Trp + 2am (^j- + a, k^ 



e [0, 



p = 2m+l {m = 0,l, 



(50) 
(51) 



The domains of the parameter k in (|50p and (|5ip are given by (|T7)) . whereas ckc G [0, K [k]). The boundaries of 
the stability regions ac — Uc (fc) are determined (see Appendix B) by the solutions to the functional equation 



L 



k sn — + ttc. A: sn ac. A: cn — + etc, fc cn ac,k 



L 



L 



L 



E [ - + ac, k ] + E [ — ~ ac, k 



L 



L 



L 



L 



L 



+sn y-j- — ac, k j cn — ac, k j dn + ac. A; ) + sn ( — + ac. A; ) cn ( — + ac. A; ) dn ( Uc, A; ) = (52) 

in the case (1501) . and to the hmctional equation 



l-fc2 



sn 



ac,k sn 



Qfc, A: I cn I — + ac, k \ cx\ 



L 



ar- k 



X { E \ J + ac,k 



L 



E[ - -ac,k] - k^ 



{^+ac,k) cn(|+a„fc) ^ sn(^ 



a, 



,k) 



cn (j - ac, k) 



dii{^ + ac,k) 



dn(^ - ac,k) 



. - ac, k) cn (-1 - ac, k) sn + ac. A:) cn + ac. A:) 



dn 



,fc) 



dn (^ - ac 



,k) 



= (53) 



in the case (|511) . The relevant solutions to ((52| and (|53|) must satisfy the conditions ac (fcp) = (p = 1,2,...). 
Making the substitution 

ac^K (k) - 7c (54) 

and proceeding to the limit A; = 1 in Eq. (f52|) with p = 0, we arrive at Eq. (|36|) . Accordingly, for A: ^ 1, the 
asymptotics of ac (A;) coincides with that of Pc {k) [relation ([57]) ]: 



7c (i), 



(55) 



where the dependence 7c = 7c (i) is represented by the graph in Fig. [5l The feature (I55p is clearly reproduced in 
Fig. m where we present the stability regions of (I50p and (j51|l obtained by numerical evaluation of (I52|) and (j53|) . 
respectively, for L = 0.3, 1,3. Moreover, as could by expected from the general arguments at the beginning of this 
section, the limiting form (A; = 1) of the current-carrying solution (jSOp for p = coincides with the limiting form 
(fc ~ 1) of the current-carrying solution ([55]) . i.e., 

lini (j)o = lim 0s 

where (pi = 4>i (y) is given by ([55)1 . 

It is interesting to note that equations ([5^ and (|53p have exact analytical solutions at the points k — k* {p — 
0, 1,2,.. .), where k* arc implicitly determined by the equations 



1 



p+- fc;/v(fc;)=L, p = 0,l,2,.... 



(56) 



Namely, 



ac (k*) 



p = 0,l,2. 



(57) 
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The role of these solutions is discussed in Sec. IV. 

Upon the substitution of ([50]) and ((5T|) into ([5]) and ([S]), we obtain, respectively. 



H = 



J = 



2k 



dn-i 



dn" 



dn" 



k - dn 



— a,k 

— a,k 
a g [0, ad 



(58) 
(59) 



for the case (f50| {p — 2m, m = 0, 1, . . .), and 

H -- 
J 



1 








dn 




2k 


J 


1 








dn 









dn 



+ a,k ] — dn 



— a, fc 

■ — a,k 
a £ [0, aj , 



(60) 
(61) 



for the case (fSTj) (p = 2m + 1, to = 0, 1, . . .), with the domains of the parameter k being determined by (|46p . (|47|) . In 
the hmit fc = 1, equations ([581), dSS]) for p = take the form ([41]), (jiS]), as they should. 

By setting a = ac 'm ([551) - (pl]) . we can obtain relevant parts of the dependence Jc = Jc (^^) for arbitrary L e (0, oo): 
this is the subject of the next section. However, we want to conclude this section by demonstrating how the above 
exact analytical results reproduce the well-knownii^ii^iii " Fraunhofer pattern" of = Jc {H) in the limiting case 



L^f «1. 



In the case L <C 1, the solutions to Eqs. (|T7)) are 

2L 

/cpw — <1, p=l,2,... 
(see Fig. m for L — 0.3). Accordingly, the domains of the parameter k [relations become 



p = : € 



2L 



1 ; p=l,2,. 



fc e 



2L 



2L 



(p + 1) TT pvr 



(62) 



(63) 



Therefore, we focus our attention on the case fc <C 1. In this limit, for k ^ kp {p — 1,2, . . .), the solution to both Eq. 
([52| and Eq. ^ is 



TT 



For fc <C 1, equations (|58l) and ([60|l yield 

whereas Eqs. ([5^ . (pT|) become 

j« (_l)P:^sin 
^ ' 2 



2L\ 



sin (2a) , a G 



0, ■ 



P = 0,l,2,, 



Combining relations (|63p - (|66p . we arrive at an approximate dependence J J {H, a) for L 



< 1: 



As a result. 



J« — |sin(i7W^)|sin(2a), a£ 
2H 



Jc{H)^—\smiHW)\. 



TT 

«'4 



(64) 



(65) 



(66) 



(67) 



(68) 



Our derivation clearly reveals limitations of the approximate relation (the "Fraunhofer pattern"): strictly 
speaking, in the field range < H < 1 (i.e., when fc — > 1), it can be regarded, at most, as a reasonable interpolation. 



12 



Moreover, the approximation ([68)) breaks down near the boundaries of the stabihty regions H x p = 1,2,..., (i.e., 
when k « kp). Unfortunately, these hmitations are not accounted for in elementary derivations^i^!^ of 

Finally, we note that, as is clear from ([65)1 . the actual expansion parameter in relations ([66)) - ([68|) is H^^ <C 1 
rather than L ^ 1. Therefore, for H ^ 1, the approximation is valid for arbitrary W S (0, cxd). (This fact 
was first pointed out in Ref.-i^.) For reference purposes, we present the corresponding (H 3> 1) asymptotics of the 
current-carrying solutions <j)p- 



(y) « pTT + 2Hy + 2a 



(-If 



2i72 47J2 



[sin {2Hy + 2a) - 2Hycos (HW) cos (2a) 



a e 



Or 



p = 0:He 1, 



W 



p= 1,2,. 



He 



pn (p + 1) TT 
W' W 



(69) 



IV. MAJOR PHYSICAL RESULTS 



A. Stability regions in the plane {H, J) and the dependence Jc — Jc (H) 

Relations ([5^ . (|^ and ([55)) - ([CT|) map the stability regions of the current-carrying solutions ([55)1 and ([50]) . ([?T|) . 
respectively, from the parametric planes (/3, fc), (a, k) to the physical plane {H, J). In Fig. [51 we present the results 
of this mapping for the data of Fig. [D As already noted (Sec. Ill), the boundaries of the stability regions represent 
the dependence Jc ~ Jc (H) that consists of an infinite number of separate branches. 

As can be easily seen, the structure of the stability regions [including the boundaries Jc = Jc (H)] is qualitatively 
the same for all the considered cases: L — 0.3 (a "small" junction), L = 1 (a "medium" junction), and L = 3 (a 
"large" junctions). Thus, as the solutions 05 [Eq. (|33p ] and 0o [Eq. ([50)) ] constitute two different branches of the 
same current-carrying solution, their stability regions (labeled by the indices s and p = 0, respectively) merge to form 
a unified stability domain. The transformation (j)s < — > 0o occurs on the internal boundary ([4T]) . ([42[) (represented by 
the dashed line in Fig. [5]), where these two solutions coincide with the elementary solution (pi [Eq. (|38p ]. 

The stability regions corresponding to 4'p=2m and 0p=2m+i interchange and form an infinite set. Significantly, for 
arbitrary L £ (0, cxd), each two consecutive stability regions, labeled by p and p + 1, overlap in the field range 

^Hj - 1 <H<Hp, (70) 

where Hp are the roots of 

{p+1)k(^±^ ^HpL, p = 0,l,2,.... (71) 

Indeed, for J — 0, the stability regions of the solutions (pp [Eqs. (|44p . ([45| ] are given by the field intervals^ 

p = 0:0<H<Ho; (72) 
P^l,2,...: ^Hl_^-1<H <Hp- (73) 

hence relation (|70p . Moreover, for sufficiently large L, the overlap may involve several consecutive stability regions: 
see Fig. [6) for L = 3. In contrast, the overlap decreases with an increase in p and a decrease in L: see the insert in 
Fig. [51 for L = 0.3. The overlap of the stability regions results in multivaluedness of the dependence Jc = Jc (H). 

For L > 1, the multivaluedness of Jc — Jc{H) was found by numerical evaluation.^ However, the fact that this 
multivaluedness is an intrinsic feature of any Josephson junction (even with L <C 1) was not noticed because of the 
absence of exact analytical solutions. 



B. Unquantized Josephson vortices 

In contrast to the case J = the discrete parameter p of the exact solutions ([501) and ([51]) for J > cannot 
be identified with the number of Josephson vortices, although relation ([48l) still holds. The reason is the occurrence 
(for certain values of H and J) of unquantized vortices carrying fractional flux $ e (i^oj^o); where $o is the flux 
quantum. (In our dimensionless units, $0 = ""•) To clarify the situation, we should consider spatial distribution of 
the local magnetic field h and of the Josephson current density j for J = Jc- 
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FIG. 6: The stability regions of 4>s and (?!>p (p = 0, 1, 2 . . .) in the physical plane {H, J) (shaded) for L — 0.3, 1, 3. The dependence 
Jc = Jc (H) is given by the solid lines. The dashed line represents the internal boundary where (ps = (po = 4>i.- 



As follows from ([T|) and jS]), the local magnetic field h obeys the linear homogeneous second-order differential 
equation 



2 -cos (/)(?/) /i. (74) 



dy 

Combining Eq. ([71]) and Eqs. PTI)) . ([TT|) for the boundary of the stability region, we obtain 



^,(L)^{L)=i;oi~L)fi-L). (75) 
dy dy 
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or 



Taking into account ([T2)) . we find that either 

The solution (t)s [Eq. ((33)) ] satisfies relation ([77]) everywhere on the critical curve Jc = Jc{H). Moreover, for this 
solution, ^ (y) > for any y € [— L, L\. [Accordingly, j {y) > for any y e [— L, L]: see © and ([3]).] 

The behavior of 0p [Eqs. ([50)1 and ([5T|l ] is more complicated. First, we note that (pp satisfy ([76)l at those values of 



i7 for which Jc = 0. This occurs bX H = Hq (for 0o) and at H = y Hp_^ — 1, i/p (for 0pwith p — 1,2, . . where Hp 
are determined by ([7T|l : as a matter of fact, this case has been considered in detail in Ref^. 

In addition, relation ([76]) is satisfied by 0p at such fields H — H* (p = 0, 1, 2, . . .) that Jc — (2i/*) ^ . These fields 
are given by 



rr* ^ 



p = 0,l,2,. 



where /c* are determined by ([56]). At H = H*, we have: (/)p (-L) = 0, (/)p (i) = tt (2p + 1). 

At the rest of the points on the critical curves Jc — Jc (H), the solutions (fip [p = 0, 1, 2, . . .) satisfy ([77]) . In particular, 
we have: (a) ^ (±L) > for H < H* (p = 0, 1, 2, . . .), because cjip {-L) g (O, |) and 0p (L) G (27rp, 27r (p + i)); (b) 
^(±L) < for iJ > iJp* (p = 0,1,2,...), because 0p G (-7r,0) and 0p (L) € (27r (p + i) , 27r (p + 1)) . 

To estabhsh the types of Josephson- vortex structures that are represented by the solutions (j)p (p = 0, 1, 2, . . .) on the 
critical curve Jc = Jc (iJ) > 0, we have to classify the points of local minima of h. Thus, for H < H* [p = 0, 1, 2, . . .), 
the first minimum is positioned at j/ = —L, where ^ {—L) > 0. The rest of the minima (for p > 0) are positioned at 
y = yn (n = 1, . . . where (j)p (?;„) = 27rn and ^ (y„) = 0. 



For H = H* (p = 0, 1, 2, . . .), the first minimum is positioned at y = —L, where (/)p {—L) = and ^ {—L) = 0. The 



rest of the minima (For p > 0) are positioned at y = y„ (n = 1, . . . where 0p (?/„) = 27rn and 4^ (t/„) = 



For iJ > Hp (p = 0, 1, 2, . . .), we have a minimum aX y = L, where ^ {L) < 0. The rest of the minima are 

positioned at y = y„ (n = 0, . . . where 4>p (yn) = 27rn and ^ (?;„) = 0. 

Bearing in mind that a Josephson vortex is located between two consecutive local minima of hj~ we arrive at the 
following physical interpretation of (f)p: 

(i) H < Hp {p = 0,1,2,...). The solution 0o represents a vortex-free configuration. The solutions labeled by 
p = 1,2, .. . represent configurations with p — 1 quantized Josephson vortices located between the points y„, yn+i 
(n — l,...,p — 1) and carrying flux $ = $0. In addition, these configurations contain a single unquantized vortex 
carrying flux $ G *I'o) and located between y = —L and y = yi; 

(ii) H = Hp [p = 0,1,2,...). The solution 0o represents a vortex-free configuration. The solutions labeled by 
p = 1, 2, . . . represent configurations with p quantized Josephson vortices located between the points y — —L, j/i, and 
Vn, Vn+i (n = 1, . . . ,p- 1; p > 1); 

(iii) H > Hp {p = 0,1,2,...). The solutions labeled by p — 1,2,... represent configurations with p quantized 
Josephson vortices located between the points y„, yn+i (n = 0,...,p— 1). In addition, all these configurations 
{p = 0, 1, 2, . . .) contain a single unquantized vortex located between yp and y = L. 

The above general analytical conclusions are illustrated in Figs. [7] and [8] For simplicity, in Fig. [71 we restrict 
ourselves to the first two critical curves Jc = Jc (H) of the junction with L — I. Spatial distribution of h and j 
at typical points 0-7 on these curves is presented in Fig. [H were we also mark the locations of both quantized and 
unquantized Josephson vortices. 

In conclusion, we want to emphasize that, as follows from continuity arguments, unquantized Josephson vortices 
persist in certain two-dimensional domains on the plane {H, J), where J < Jc. Therefore, the existence of such 
vortices is a typical feature of any Josephson junction in the presence of externally applied magnetic fields and 
transport currents. 
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FIG. 7: The first two critical curves Jc = Jc {H) for L = 1. Spatial distribution of h and j (Fig. |8]) is evaluated at typical 
points 0-7: (0) H = 0.00, Jc = 0.99; (1) H = 0.43, Jc = 0.87; (2) H = 0.61, Jc = 0.76 (tps = 00 = 0;.); (3) H ^ Ho = 0.81, 
Jc = (2//o)~^ = 0.61; (4) H = 1.60, Jc = 0.05 (the first curve); (5) H = 1.60, Jc = 0.05 (the second curve); (6) H = HI = 2.36, 
Jc = {2Hl)-^ = 0.21; (7) Jf = 3.01, Jc = 0.05. 



C. Generalizations 



The restriction > 0, J > imposed at the beginning of Sec. Ill can be easily removed. Physical solutions that do 
not obey these restriction are expressed via the solutions 0^, (/)p and (pi [Eqs. ([55]) . ([SO)) . ([5T|) . and respectively] 
by means of elementary symmetry relations. 

1. The case H <0, J>0: 



2. The case H>0, J <0: 



3. The case H <0. J <0: 



,/3->0s,-/3; (j)p,a -(t>p,-a; 0;, 7 -> 27r - 0;, -7. 



,/3 ^ -</>s, -/3; '/'p,a ^ 0p,-a; 7 ^ 0; - 27r, -7. 



(78) 
(79) 
(80) 



Finally, we note that the consideration of this paper equally applies to a generalized form of the boundary conditions 
that takes into account possible asymmetry in the injection of the transport current, namely)^i^ 



dy 



(±L) = 2{H±a^J), 



(81) 



or, equivalently. 



^^2 



J 



1 



dy dy 



where 

a± > 0, a_ + a+ = 1. 
The effect of the generalized boundary conditions (|8ip is illustrated in Fig. [SI 
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FIG. 8: Spatial distribution of h (solid line) and j (dashed line) for points 0-7 in Fig. [7] The location of Josephson vortices 
is marked by vertical arrows: the dashed arrows correspond to unquantized vortices [figures (4), (5) and (7)]; the solid arrows 
correspond to quantized vortices [figures (6) and (7)]. 
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V. SUMMARY AND CONCLUSIONS 



Summarizing, we have derived the complete set of exact physical solutions to the general boundary-value problem 
©, (El): 05 [Eq. (O] and (l)p {p = 0, 1,2, . . .) [Eqs. dSO]), dST])] complemented by the symmetry relations ([78l) -(l80l). 
The obtained solutions describe the current-carrying states of the Josephson junction of arbitrary length W = 2L E 
(0, oo) in the presence of an externally applied magnetic field H € [—He, He), where He is the thermodynamic critical 
field of the superconducting electrodes. The most direct application of these solutions is straightforward evaluation of 
the dependence Jc = Jc (H) (for arbitrary W and an arbitrary mode of the injection of J) by means of the algorithm 
of Sees. Ill and IV. 

Mathematically, 4>s and 4>p (p = 0, 1, 2, . . .) represent the complete set of particular solutions to ^ that are stable 
under the condition that ^ is fixed at the boundaries y — zLL, and they possess a number of interesting properties. 
For example, the solutions 0s and 0o constitute two different branches of the same stable solution: for fc — ^ 1, both of 
them turn into the elementary solution (pi [Eq. ([55]) ]. In physical literaturefi^i^ the elementary solution 0; is usually 
identified with "the Josephson vortex". Indeed, if Eq. ([3]) were considered on the infinite interval (— cx),cx)), this 
solution would be nothing but the well-known^"^ static soliton of the sine-Gordon equation, positioned aX y — j and 
stable for arbitrary I7I < 00. However, on the physically realistic finite interval [— L, L], the solution (jji proves to be 
stable only for I7I > 7c, where jc is determined by Eq. ()36|) . and, as shown in Sec. IV, it has nothing to do with any 
vortex (or soliton) configurations. 

As could be anticipated, the exact analytical solution of the problem that remained unresolved for over four decades 
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has revealed some unexpected physical features. For example, contrary to a wide-spread beliefjii^i^ it clearly shows 
that there is no qualitative difference between Josephson junctions with W ^ 1 and those with W 1. Thus, the 
exact analytical dependence Jc = Jc (H) proves to be multivalued even for arbitrarily small W. Therefore, hysteresis 
is an intrinsic feature of any Josephson junction with W € (0, cxd). 

However, we think that the most important physical conclusion that can be drawn from the exact analytical solution 
is the existence of unquantized Josephson vortices. Indeed, recently, the possibility of finding unquantized vortices in 
different types of superconducting systems (including Josephson ones) has attracted considerable attention: see, e.g., 
Refs.^^i^ and references therein. In most theoretical models, unquantized vortices appear as a result of unconventional 
properties of the superconductors themselves, such as, e.g., the existence of two superconducting order parameters;^ 
d-wave pairing combined with the inhomogeneity of grain boundarieSf^^ etc. By contrast, we have shown that the 
presence of unquantized Josephson vortices near the external boundaries is a typical feature of any classical Josephson 
junction, provided the transport current J is sufficiently close to Jc for certain finite values of H . 

From a mathematical point of view, it would be desirable to know whether the quantity Jmax discussed in the 
Introduction indeed coincides with Jc evaluated in this paper. Although we have been unable to find a general 
analytical proof, our detailed comparisons with the numerical results of Refs.-'U, (not presented here for brevity 
reasons) imply that the identity Jmax = Jc can be accepted without reservation. 

Finally, we want to remind once again that Eq. ([3]) is just the static version of the well-known sine-Gordon 
equation [Eq. (jAl[) with k = 0]. Given that the sine-Gordon equation finds a lot of applications in condensed-matter 
and elementary-particle physics^ we hope that our exact analytical solution may find applications outside the field 
of superconductivity as well. 
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APPENDIX A: ALTERNATIVE FORMULATION OF THE STABILITY PROBLEM 

The stability of a given solution = 0(?/) to ([3]), (jl]) can also be analyzed using the general time-dependent 



equation?. 



d^<j) ^ d(j) d'^cj) . , ^ 

^ + 2.^-^+-^-0, (Al) 
where k > 0, t > and y G {—L, L), under the boundary conditions 

^ (t, ±L) = 2H±J. (A2) 
According to linear stability theoryfiS we should seek solutions to (|Aip . (|A2p in the form 

4>{t,y) = cbiy) + e--*Ciy), (A3) 

where 

max|C(y)| < 1, 

V 

and 

dc 



(±i)=0. (A4) 
dy 



Substituting ()A3|) into ()A1|) and dropping nonlinear terms, we obtain 



dy 

Equation (jASP under boundary conditions (jA4D immediately yields 



2 cos (/) (y) C = cr (2k - cr) C (A5) 



tT„± = K ± - n = 0,1,2,..., (A6) 
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where < A^i < A'2 < ■ ■ • are the eigenvalues of the Sturm-Liouville problem ([8]), ((91). Thus, we arrive at the following 
classification of stability properties of the solution (f> = (j>{y): 

i) fiQ > 0, Re (7n± > (n = 0, 1, 2, . . .): the solution is exponentially stable; 

ii) /io < 0, (To- = K — ^y^?+\Jlo\ < 0: the solution is unstable; 

iii) fiQ = 0, (Tq- = 0: the solution is at the boundary of the stability region. 



APPENDIX B: SOLUTION OF THE LINEAR BOUNDARY- VALUE PROBLEM FOR Ato = 0, i>o = V'o 

To solve the linear boundary- value problem (fTUl) - P^ . we should first find the general solution to (|10p . ft can be 
written in the form 

^0 (y) - Cixi (y) + C2X2 iy) , (Bi) 

where xi, X2 are linearly independent solutions to (fTO)) . and Ci, C2 are arbitrary constants. As to xi, we can choose^ 

XI = Cf, (B2) 
dy 



where C is a normalization constant. The linearly independent solution X2 is determined by the well-known relatioi 



X2 = Xi/^. (B3) 
J XI 

In the simplest situations, we have either 

xi {-y) = xi (y) , X2 {-y) = -X2 {y) , (B4) 

or 

xi {-y) ^ ~~xi (y) , X2 {-y) ^ X2 (y) ■ (B5) 

In the case (jB4[) . which corresponds to (p = (f>p [Eqs. ((44|) . (|45)) ] with k = kp, we have C2 = 0, i/io = Xi- Conditions 
PTjl result in Eqs. (gT]).*^ Condition ([T2]) is fulfilled automatically. 
In the case (jB5[) . which corresponds to = 0^ [Eqs. (HJ)] with /c = we have Ci = 0, 

^0 (y) = X2 (y) ^ (y^ fcc) + (1 - y] - cn {y, fe,) . (B6) 

The substitution of (|B6P into (fTT|) yields Eq. (I24p . Condition (fT2|) singles out the solution presented in Fig. [21 

Consider now the general situation, when the fimctions xi, X2 do not obey either (|B4[) or (jBSp . and, accordingly, 
Ci 7^ 0, C2 7^ 0. Upon the substitution of (|Bip into (fTTI) . we get a system of algebraic equations for Ci, C2, 

Ci^(L) + C2^(L) = 0, 
ay dy 

dXi , r^ , ^ '^X: 

with 



C.^(-X) + C2^(-L)=0, 
dy 



dy dy dy dy 



being the solvability condition. Equation (IB7p . under condition ()12p . determines the sought solution. 
In particular, in the case of (|33p with /3 = /3c, we have 

sn(?/ + /3c,fc) 
^^^^^= dn(, + /3.,fe) ' 

X2 (y) = dn(U£fc) ^"^ ^ " " ^ ■ ^^^^ 
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The substitution of (|B8| . (|B9|l into jBT)) leads to Eq. ((341). Condition ((T2]) leads to the boundary condition /?c (fee) = 
for the domain (|23p . where fc^ is determined by Eq. (|24p . 

In the case of (f50| with a — ac, the functions xi, X2 are given by 



Xi (y) = dn 

X2 (y) dn"^ + Qfc, fc) S + Qfc, fc) . 

Substituting (|B10|) and (jBlip into ([13), we get Eq. Analogously, for ^ with a = ac, 

Xi (y) =dn(^|+ac,fc) 



X2 (y) 



dn(| + Q;e,fc) 
l-fc2 



E 



(|+a„fc)- 



fc^sn (I + ac, k) cn (| + ac, k) 



dnO 



(BIO) 
(Bll) 

(B12) 
(B13) 



with Eq. (|53p being the result of substitution into (|B7p . Condition (|12p leads to the boundary conditions ac {kp) = 
(p = 1, 2, . . .) for the domains where fcp are determined by Eqs. ([37]). 
Finally, in the case of the elementary solution (|38p with 7 = 7c, 



Xi (y) = cosh ^ (y - 7c) , 



X2 (y) 



1 



2 cosh (y - 7c) 



sinh (y - 7c 



y 



(B14) 
(B15) 



Substituting (jB14p . (|B15p into (|B7p . we arrive at Eq. ((35)) that has been obtained in the main text as the limiting 
form of Eqs. §^ and dH]) (with p = 0) for ^ 1. 



APPENDIX C: EXPLICIT EVALUATION OF = FOR H > 1 

For the lowest eigenvalue fi = fiQ, the Sturm-Liouville problem ([5]), (O becomes 



+ cos(/)(y)-0o = Mo^o, y<^{-L,L), (CI) 
V'o(y)^0, ye[-L,L]. (C3) 



In the general case, the solution to (|Cip - (|C3p can be obtained using the fact that Eq. (jCip is reducible to Lame's 
equation.— However, for H ^ 1, the eigenvalue /io can be explicitly evaluated by elementary methods. 
We will seek the solution to (|Cip - (jC3P in the form of asymptotic expansions 

V'o(y)«E^o"^(y)' mo«Ea^o"^ (C4) 

n>0 n>l 

where "^q"^ and /Xq"-* are of order . Besides, we will employ the exact integral relation 

/f ^ dyVo (y) cos (y) 

Mo = X (C5) 

/_!, dyipQ (y) 

that immediately follows from (jCip - (jC3p . 

Introducing a new variable u = Hy, we rewrite (jCll) - (|C3p as 

^ + 77? [cos </-(") -Mo] V^o = 0, ue{-HL,HL), (C6) 



du2 i/2 



^,„.).0, ,C7) 

V'o (w) ^0, u G [-iJL, iJL] (C8) 
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and note that \fio\ < 1.^ Thus, the problem for ipQ ' has the form 




(Cll) 



(CIO) 



(C9) 



The solution to ((C9l) - ((Cn|) is 



ipQ {y) = const. 



(C12) 



Using relation (|C5|) . solution (jC12p and the asymptotic expansions for (pp [relation ((69|) ]. we find 



fJ-o {H, a) « {H, a) 



(-ly 



sin {HW) cos (2a) . 



(C13) 



HW 



As can be easily seen, expression (|C13p stands in full agreement with the general results of Sec. III. 
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